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Abstract 

We compare and contrast the long-time dynamical properties of two individual-based models 
of biological coevolution. Selection occurs via multispecies, stochastic population dynamics with 
reproduction probabilities that depend nonlinearly on the population densities of all species resident 
in the community. New species are introduced through mutation. Both models are amenable to 
exact linear stability analysis, and we compare the analytic results with large-scale kinetic Monte 
Carlo simulations, obtaining the population size as a function of an average interspecies interaction 
strength. Over time, the models self-optimize through mutation and selection to approximately 
maximize a community fitness function, subject only to constraints internal to the particular model. 
If the interspecies interactions are randomly distributed on an interval including positive values, the 
system evolves toward self-sustaining, mutualistic communities. In contrast, for the predator-prey 
case the matrix of interactions is antisymmetric, and a nonzero population size must be sustained 
by an external resource. Time series of the diversity and population size for both models show 
approximate 1/f noise and power-law distributions for the lifetimes of communities and species. 
For the mutualistic model, these two lifetime distributions have the same exponent, while their 
exponents are different for the predator-prey model. The difference is probably due to greater 
resilience toward mass extinctions in the food-web like communities produced by the predator- 
prey model. 
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I. INTRODUCTION 



Traditionally, problems in ecology and evolution have been addressed at very different 
levels of resolution. Typically, ecological problems are studied on a timescale of generations 
and often at the level of individual organisms, while issues in evolution are considered on 
much longer, often geological, timescales and usually at the level of species or higher-level 
taxa. However, in recent years it has been reco gniz ed that processes at the ecological and 
evolutionary scales can be strongly linked jlG, 5^, 5^, G^]- Several models have therefore been 
proposed, which aim to model the complex problem of coevolution in a fitness landscape 
that changes with the composition of the community, while spanning the disparate scales 
of both temporal and taxonomic resolution. Early steps in this direction were simulations 
of parapatric and sympatric speciation [13i] and the coupled NK model with po pul ation 
dynamics More recent contributions include the Webworld model 




Tangled-nature model Jiol . ITll . 2^ and simplified versions of the latter [s^, [H, [gI] , as well 
as network models [sl, [QfT Recently, large individual-based simulations have also been per- 
formed of parapatric and sympatric speciation [3, 25 1 and of adaptive radiation {26|]. Many 
of these models are deliberately quite simple, aiming to elucidate universal features that are 
largely independent of the finer details of the ecological interactions and the evolutionary 
mechanisms. Such universal features may include lifetime distributions for species and com- 
munities, as well as other aspects of extinction statistics, statistical properties of fluctuations 
in diversity and population sizes, and the structure and dynamics of food webs that develop 
and change with time. By changing specific features of such simplified models, one hopes 
to learn which aspects influence the observed properties of the resulting communities and 
their development with time. 

In this paper we compare and contrast two stochastic coevolution models that combine 
features from some of those mentioned above. In each, selection is provided by an individual- 
based population dynamics, while new species are provided by a low rate of mutation. The 
models are studied both analytically by linear stability theory, and numerically by large- 
scale kinetic Monte Carlo simulations. The first model (Model A) allows direct mutualistic 



interspecies interactions, and some of its properties were discussed previously 50^ 53, 65|. 
It is a simplified version of the Tangled-nature model of Jensen and coworkers 10|, lUl, l29 |. 
The second model (Model B) is a predator-prey model. For both models we obtain exactly 
the fixed-point mean population sizes and stability properties for any given community of 
species in the limit of vanishing mutation rate, using linear stability theory. These results 
enable us to define a community fitness function, cubic in the mean total population size, 
that is maximized at the fixed point for a given community of species. The fact that this 
much information can be obtained analytically makes these models uniquely well suited as 
benchmarks for more elaborate, but less tractable, nonlinear models. 

The analytical results are followed by numerical simulations of the dynamics of both 
models for nonzero mutation rates. We focus on very long simulations in a regime where 
both diversity and population size are statistically stationary, albeit with very large, strongly 
correlated fluctuations and an intermittently vigorous turnover of species. We thus study the 
intrinsic dynamics of extinction and origination of species and communities in the absence 
of external perturbations. Understanding of these intrinsic fluctuations in the station ary 
state should enable one to estimate the community's response to external perturbations 5^ 



in a way analogous to the fluctuation-dissipation relations of statistical mechanics [4^ 

Three main conclusions emerge from this combined analytical and numerical investiga- 
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tion. 

• Mutations enable the models to evolve communities that tend to maximize the commu- 
nity fitness function, subject only to constraints that are internal to the specific model. 
Although there is vigorous turnover of species and communities, all communities that 
persist for a significant time remain in the vicinity of this maximum. 

• The simulated systems exhibit power-law distributions in the lifetimes of individual 
species, as well as of communities, and the fluctuations in diversity and population 
size are characterized by approximate 1// noise. 

• The two models show distinct differences in the community turnover, which are re- 
flected in differences between the power-law exponent of the community-lifetime dis- 
tribution. In Model A, extinctions of species tend to be highly synchronized, such that 
a whole community collapses in a mass extinction on a relatively short time scale. In 
Model B, on the other hand, extinctions are more often limited to a subset of the 
resident species. This can be explained by differences between the structures of the 
interaction networks characterizing long-lived communities in the two models. For 
Model B, this network amounts to a simple food web. 

The rest of this paper is organized as follows. The models are introduced in Sec. [Tll In 
Sec. mil we perform a full analytical linear stability analysis, and the results are compared 
with large-scale kinetic Monte Carlo simulations. The simulations are described in further 
detail in Sec. HV] where lifetime distributions and power spectral densities also are reported. 
A concluding summary is presented in Sec. El Some technical details of the derivation of the 
community fitness function are given in Appendix [Xj an analytical result for the optimum 
population size of Model A is obtained in Appendix [Bl and a discussion of the consequences 
of changes in several of the model parameters is found in Appendix O 



II. THE MODELS 



In both models, selection is provided by the reproduction rates in an individual-based, 
simplified, stochastic population-dynamics model with nonoverlapping generations. This 
interacting birth-death process is augmented to enable evolution of new species by a mutation 
mechanism. The mutations act on a haploid, binary "genome" of length L, as introduced 
in Eigen's model for molecular evolution [l^, |23]- This bit string defines the species, which 
are identified by the integer label / G [0, 2-^ — 1]. Typically, only a few of these 2^ potential 
species are resident in the community at any one time. 

Individual organisms of species / reproduce asexually at the end of each generation, each 
giving rise to F offspring individuals with probability P/ before dying. With probability 
(1 — P/), they die without offspring. No individual thus survives beyond one generation. (For 
simplicity, the fecundity F is assumed fixed, independent of both species and individual.) 

During reproduction, each gene in an offspring individual's genome may undergo mutation 
(0 — i> 1 or 1 ^ 0) with a small probability, /i/L. The mutation thus corresponds to 
diffusional moves from corner to corner along the edges of an L-dimensional hypercube 



|22l . |23| |. A mutated individual is assumed to belong to a different species than its parent, 
with different properties. Genotype and phenotype are thus in one-to-one correspondence in 
these models. This is clearly a highly idealized picture, and it is introduced to maximize the 
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pool of different species available within the computational resources. The approximation 
is justified by a large-scale computational study of a version of Model A, in which species 
that differ by as many as L/2 bits have correlated properties [13] • Quite remarkably, this 
study reveals that the correlated model has long-time dynamical properties very similar to 
the uncorrelated Model A studied here. 

The reproduction probability Pi{t) for an individual of species / in generation t depends 
on the individual's ability to utilize the amount of external resources available, i?, and on 
its interactions with the population sizes n,j{t) of all the species present in the community 
at that time. The dependence of Pi on the set of nj is determined by an interaction matrix 
M js^ with elements Mjj G [—1,1] in a way defined specifically in the next paragraph. 
We emphasize that M is chosen randomly at the beginning of each simulation run and is 
subsequently kept constant throughout the run. If Mjj is positive and Mj/ is negative, then 
/ is a predator and J its prey, and vice versa. If both matrix elements are positive, the 
species interact directly in a mutualistic way, while both elements negative implies direct 
competition. 

Specifically, the reproduction probability for species J, Piit), depends on R and the set 
{nj{t)} through the nonlinear form, 

^'^^^ = l + exp[-Aj{R,{nj{t)})] ' 
with the density dependent argument 

Ai{R, {nj{t)}) = -hi + ViR/Ntot{t) + J2 Mijnj{t)/Ntot{t) - Ntot{t)/No . (2) 

J 

Here bj can be seen as the "cost" of reproduction (always positive), and r]i (positive for 
primary producers or autotrophs, and zero for consumers or heterotrophs) is the ability of 
individuals of species / to utilize the external resource R. The latter is renewed at the same 
level each generation. It does not have independent dynamics. The total population size is 
Ntot{t) = ^j^j(t), and the constant A'o is an environmental carrying capacity that 



prevents Ntotit) from diverging to infinity. It may be seen as representing implicit resource 
limitations not explicitly included in R, such as available space. 

For large positive Aj, the individual almost certainly reproduces, giving rise to F off- 
spring, while for large negative A/, it almost certainly dies without offspring. The repro- 
duction probability Pj, together with its argument A/, play the role of a functional response 



for this class of models |16l . |37|. The two externally determined parameters that influence 
the population size, the resource term R and the carrying capacity A'^o, play very different 
roles. From Eq. ([2]) it is seen that R encourages population growth, especially for small total 
population sizes, while it has very httle effect for large Atot- It can thus be thought of as 
representing a finite amount of available food. The carrying capacity A"o, on the other hand, 
always discourages growth, but its effect is only appreciable for large population sizes. It 
cannot maintain a nonzero population by itself and is thus best thought of as representing 
an overall limitation, such as finite available space. If i? = 0, a nonzero population can 
therefore only be maintained by mutually positive interspecies interactions. This rather 



unrealistic aspect is shared by the Tangled-nature model [10|, lUl, 129 . 

In this work, the model parameters are chosen to represent the realistic situation that the 
number of species resident in the community at any time is much smaller than the number 
of potential species (i.e., that Af{t) < 2^), and also that M{t) < A'tot(t). 
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An analytic approximation describing the time development of the mean population sizes 
(averaged over independent realizations), {ni{t)), can be written as a set of coupled difference 
equations, 

Mt+1)) = {nj{t))FPj{R,{{nj{t))})[l-fi] 

+ {f,/L)Fj2{nKiimPKii){Rd{njm) + 0{fi') , (3) 

K{I) 

where K{I) is the set of species that can be generated from species / by a single mutation 
("nearest neighbors" of / in genotype space). 

A. Model A 



Model A was introduced and studied in [50|, l65[. In this model, the Mjj for I ^ J are 
stochastically independent and uniformly distributed on [— 1,+1], while the intra-species 
interactions Mu = 0. The external resource R and the reproduction costs bj are equal to 
zero, and the total population size Atot(^) is limited only by the carrying capacity Nq. The 
model is found to evolve through a succession of quasi-stable, mutualistic communities. 

B. Model B 

Model B is a predator-prey model. This is implemented by making the off-diagonal 
part of M antisymmetric. In order to keep the connectance of the resulting communities 
consistent with food webs observed in nature [H, ill, the {Mij,Mji) pairs are chosen 



nonzero with probability c = 0.1. The nonzero elements in the upper triangle of M are 
chosen independently and uniformly on [—1, +1]. This model does not include a population- 
limiting carrying capacity [i.e., formally, Ai'o = oo in Eq. ([2])], and the community is supported 
by a constant external resource, R. Only a proportion p of the 2^ potential species are 
producers that can directly utilize the resource (for the numerical data reported here, we 
use p = 0.05). Thus, with probability (1 — p) the resource coupling rji = 0, representing 
consumers, while with probability p the rij are independently and uniformly distributed on 
(0, +1], representing producers of varying efficiency. In addition to these constraints on M, 
we require that producers {rjj > 0) always are the prey of consumers {rjj = 0). In other 
words: the case 777 > and f]j = with Mjj = —Mjj > is forbidden. Whenever it 
occurs, this situation is corrected during the construction of M by reversing the signs of 
Mjj and Mjj for the pair in question. Consumers at higher trophic levels are allowed, and 
they are indeed observed in the resulting communities (see further discussion in Sec. II Vp . 
The population sizes are limited by independent reproduction costs bj that are uniformly 
distributed on (0, +1], and by negative intra-species interactions Mjj independently and 
uniformly distributed on [—1,0). The model is found to evolve through a succession of 
quasi-stable, food-web like communities. Some preliminary numerical results were presented 
in ET 
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III. LINEAR STABILITY ANALYSIS 



The form of A/ in Eq. ([2]) represents frequency- dependent interactions that describe 
universal competition and absence of adaptive foraging, and so is not very reahstic. However, 
the normahzation by A^tot (t) has the advantage that it turns Eq. with fi = into a set of 
hnear equations for the average population sizes in the fixed-point community corresponding 
to a particular set of species. The set of equations can be solved analytically to give the 
average total population size, as well as the average population size of each species and 
the stability properties of the community. This mean-field level analysis also enables us 
to construct a community fitness function that for a given set of species is maximized by 
the fixed-point community, and whose width is proportional to the size of the population 
fluctuations around the fixed point, caused by the statistical birth-death process. (For a 
detailed discussion of the fluctuations in Model A in the absence of mutations, see jGSj.) 
For small mutation rates, this picture remains valid as a description of the population 
fluctuations on relatively short time scales, where the species composition of the community 
remains constant, except for small populations of unsuccessful mutants. 



A. Fixed-point communities 

To obtain the stationary solution of Eq. ([2]) with /i = for a community of Af species, we 
require Pj = 1/F for all Af species. Equations ([T]) and ([2]) then yield the Af linear relations 

where bj = bj — \n{F — 1). (For simplicity, we have dropped the ( ) notation for the 
average population sizes. The asterisk superscripts denote fixed-point solutions.) In a vector 
notation where w is a column vector and its conjugate row vector, n* is the column vector 
composed of the Af nonzero n*j, while 1^ is an A/"- dimensional row vector composed entirely 
of ones. Thus, the total population size is given by N^^^ = ^ju} = P^n*, and Eq. (jl]) takes 
the matrix form 

- 6iV,;, + ffR + Mn* - r(iV4J ViVo = . (5) 

Here, b is the column vector whose elements are bj, ff is the column vector whose elements 
are rji (in both cases including only those Af species that have nonzero n*j), M is the 
corresponding Af x A/ submatrix of M, and 1 is an A/"- dimensional column vector of ones. 
The solution for n* is 

(6) 

where is the inverse of M. (See below for a discussion of the effects of a singular M.) 
To find each n*j, we must first obtain iV^^^ = Vn*. Multiplying Eq. from the left by V, 
we obtain the quadratic equation for A^tot, 

R£ + eiv;, - (iv;j ViVo = . (7) 



n 



-M 



fJR - bN,l, - i(iv,;j ViVo 
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FIG. 1: (Color online.) Solutions for (a) the total fixed-point population size N^^^ and (b) the 
maximum community fitness function <I>max = '^'(-^tot)- Both are shown vs Q for fixed £ = 0.61. 
The solid vertical lines represent the absolute maximum value for Q in Model A, while the dotted 
vertical lines represent the analogous limit for Model B. The points marked "Model A" and "Model 
B" were obtained from large-scale Monte Carlo simulations as described in the text. The calculation 
of the point marked "Model A theoretical" is described in Appendix iBl 



The coefficients, 

e^lZ^ a„d e^Z^L, (8) 

can be viewed as an effective interaction strength and an effective coupling to the external 
resource, respectively. Approximate expressions for B and S that are less accurate but more 
intuitive are obtained in Appendix |Al The nonnegative solution of Eq. ([7j) is 



^t;t = ^ + y J +B£No. (9) 

Figure [11(a) shows N^^^ as a function of B for two choices of Nq and R at fixed S. Special 
cases of the solution are 

{BA^o for = and B > 
for i? = and B < 

y/R£I% for B = and ^ > " 
-RS/Q for A'o = oo and/or l^M-^l = 

To find each n*j separately, we now only need to insert A^^^ in Eq. ([6]). 

Only those n* that have all positive elements can represent a feasible community 5l| . If 
M = or is otherwise singular, the set of equations ([5]) is inconsistent for A/" > 1, unless 
bj and rjj both are independent of /. The only possible stationary community then consists 
of one single species, the one with the largest value of —bj for Model A or the one with the 
largest rjj/bi for Model B. This is a trivial example of competitive exclusion 0, H, [s^], and 
stable multispecies coexistence in these models requires a nonsingular interaction matrix 
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FIG. 2: (Color online.) The community fitness function shown vs A'tot for Model A (black) 
and Model B (gray, red online). For Model A the parameters are: F = 4:, R = 0, and Nq = 2000. 
For Model B they are: F = 2, R = 2000, £ = 0.61, and Nq = oo. The values of 9 used are 1.61, 
corresponding to the average value taken by Model A in the simulations, and —0.15, corresponding 
to the average value taken by Model B. <I>(A'tot) has a nontrivial maximum for Model A at positive 
(solid black curve marked A), and for Model B at negative G (solid gray curve marked B (red 
online)). The circular and square data points are the results of Monte Carlo simulations and 
correspond to the equally shaped points in Fig. [H They lie very close to the maximum of <I>(A'^tot) 
for each model. The dashed curves show the physically inaccessible cases of Model A with G < 
(stable, absorbing state at A'^tot = 0) and Model B with G > (monotonically increasing $(A'tot))- 



Equation ([7]) can be seen as a maximization condition for a community fitness function, 



3No 



(11) 



whicfi is cubic in A^tot 67|]. Tfie dependence of $ on A^tot is shown in Fig. [2] for Models A and 
B at two different values of B. In each model, the maximum of $ represents the fixed-point 
value of iVtot, while its width gives a measure of the extent of the fluctuations about the 
fixed-point population in a community of fixed composition [65|. The maximum value of $ 
with respect to A^tot, '^'max = "^(Atot); is shown in Fig. [T](b) vs 9 for two different values of 
A'o at fixed S. 

From Eq. (ITOl) and Fig. [T](a) it is seen that Model A (finite A^o and R = 0) undergoes 



12, 57 



68 



at e = 0. For 6 < 0, A^tot 
Such behavior function of 



a transcritical bifurcation or exchange of stabilities 
vanishes, while for B > it increases linearly with 
a control parameter is common in many nonequilibrium systems. In addition to population 
dynamics and the logistic map (to which the present models withu = are closely related), 
these also include lasers and autocatalytic chemical reactions [28|, l57| . On the other hand, 
for Model B (A'o = oo and R> 0), A^tot diverges to infinity as 9 approaches zero from below, 
and it would be infinite for > 0. 

The results discussed in the previous paragraphs would not be of evolutionary relevance if 
B were just an externally fixed control parameter, as it would be in the absence of extinctions 
and mutations. However, mutations enable the community not only to maximize $(A"tot) for 
fixed parameters, but also to change the parameters in response to new mutations, leading 
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to the possibility of further increasing $max- Numerical results from the large-scale kinetic 
Monte Carlo simulations discussed in Sec. II VI show that, on average, the communities initially 
progress toward, and later settle down to fluctuate near, the maximum value of 6 (and 

thus the maximum value of $max) compatible with the constraints on M and b. These 
constraints depend on the speciflc model as follows. For Model A, which allows direct 
mutualistic interactions (i.e., Mjj and Mjj both positive), can be positive, thus enabling 
communities with nonzero population size, even for R = 0. For the predator-prey Model B, 
on the other hand, the antisymmetric structure of M forces B to be nonpositive. A detailed 
discussion of the time evolution of the mean-field parameters, Q, S, N^^^, and $max5 is given 
for both models in Sec. IIV Al 

In Monte Carlo simulations of Model A (in which the off-diagonal elements of M are 
uncorrelated and uniformly distributed on [—1,-1-1]), it was found that, after the initial 
period characterized by an average positive trend in the mean-field parameters, the com- 
munity spent most of its time in a succession of quasi-steady states (QSSs), separated by 
brief bursts of intense evolutionary activity jsO]. AH the QSSs studied in detail were found 
to be mutualistic, with M/j = 0.78 ± 0.03. [Here, the overbar represents averages over all 
the ten QSSs listed in Table I of [50|.] The average of A^tot, taken over all the 16 realizations 
of 2^^ generations that were studied, was Atot = 3201 ± 8. (The average over only the ten 
QSSs agrees with the total average to within the statistical errors, showing that the periods 
when the system is not in a QSS contribute negligibly to the overall time averages.) From 
these averages we can use Eqs. ffTOl) and f|TT]) to estimate the average mean-field quantities, 
e = 1.61 ± 0.01 and = (2.09 ± .06) x 10^ In Fig. W^a), A^ is shown vs 9 as a black 
dot, while the corresponding value of $max is shown the same way in Fig. [It^b). 

The results for Model B that are included in Fig. [T] were also obtained from Monte Carlo 
simulations of 2^^ generations. (See Sec. IIVI for further details.) The parameter values used 
in the figure were extracted as the average values from thirteen QSS communities identified 
at late times in twelve independent simulation runs: S = 0.61 ± 0.04 and G = —0.15 ± 0.01. 
The value of the total population size given in the figure is the total time average over the 
simulations, averaged over all twelve runs, Atot = 7726 ± 303. (Like for Model A, averages 
over only the particular observed QSS communities agree with this total average to within 
the error bars.) The large uncertainty in Atot is a direct consequence of the steep slope of 
A^tot versus B closely below G = for Model B (see Fig. [TJ^a)). From these estimates and 
Eq. we further get the estimate $max = (2.49 ± 0.18) x 10^. 



B. Stability of fixed-point communities 



6ij + A 



ij 



(12) 



The internal stability of an AT-species fixed-point community is obtained from the matrix 
of partial derivatives, 

dni{t + l) 
dnj{t) 

where Su is the Kronecker delta and A/j are elements of the community matrix A (39| . 
Straightforward differentiation yields 

1 



A 



IJ 



1 



N* 



Mjj- 



Rr]i + (Mn* 

N* 

-'^tot 



1 

Ao 



(13) 
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Time (generations) Time (generations) 



FIG. 3: (Color online.) Typical time series of the mean- field, fixed-point parameters for Model A 
(A^o = 2000, R = 0, F = 4, L = 13, and n = 10~^), emphasizing the transient early-time regime by 
using a logarithmic time axis, (a) The effective interaction strength 0. The solid, horizontal line is 
the absolute upper bound on Q for Model A, (1 + ln3) ~ 2.10, while the dashed, horizontal line is 
the most probable value, 1.618. Both are derived in Appendix |Bl (b) The fixed-point population 
size A^*ot and the corresponding inaximiim VclIug of th.6 comiimnity fitness function, ^max* 

The 

dashed, horizontal lines are the corresponding values obtained from the most probable value of Q. 
For comparison, the simulated total population size (sampled on a much finer time scale) is shown 
light gray (green online) in the background. 



where (Mn*)/ is the element of Mn* corresponding to species /. In order for deviations 
from the fixed point to decay monotonically in magnitude, the magnitudes of the eigenvalues 
of the matrix of partial derivatives in Eq. (fT2l) . A + 1 where 1 is the A/"- dimensional unit 
matrix, must be less than unity. The values of the fecundity F that are used in this work 
(4 for Model A and 2 for Model B) were chosen to satisfy this requirement for Af = 1. 

Since new species are created by mutations, we must also study the stability of the fixed- 
point community toward "invaders." Consider a mutant invader i. Then its multiplication 
rate, in the limit that rii <^ nj for all Af species J in the resident community, is given by 

^-(^+^) = Z (14) 

n,(t) l + exp[-A,(i?,K})] ■ ^ ^ 

The Lyapunov exponent, ln[nj(t -|- l)/ni(t)], is the invasion fitness of the mutant with 



respect to the resident community [ij, l38|. A characteristic feature of the QSS communities 
observed in both Model A and Model B is that only about 1% of the mutants that are 
separated from the resident community by a single mutation ("nearest-neighbor species") 
have multiplication rates above unity (and most of those are between 1.0 and 1.1). In fact, 
this is true only to a slightly lesser degree for mutants separated by two or three mutations 
from the resident species ("next-nearest neighbors" and "third-nearest neighbors") (50| . 
Thus, a string of rather unsuccessful, or at best neutral, mutations is necessary to bring 
significant change to a QSS community - a fact that to a large extent accounts for their 
high degree of stability. 
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Time (generations) Time (generations) 



FIG. 4: (Color online.) Typical time series of the mean-field, fixed-point parameters for Model B 
(A'o = oo, R = 2000, F = 2, L = 13, and /i = 10~^), emphasizing the transient early-time regime 
by using a logarithmic time axis, (a) The effective interaction strength and the effective coupling 
to the external resource, £. The solid, horizontal line at zero is the absolute upper bound on O 
for Model B, while the dashed, horizontal lines are the averages over late-time communities, also 
shown as points in Figs. [T] and O (b) The fixed-point population size N^o^ and the corresponding 
maximum value of the community fitness function, 'I'max- The dashed, horizontal lines are obtained 
the same way as in (a). For comparison, the simulated total population size (sampled on a much 
finer time scale) is shown light gray (green online) in the background. 

IV. COMPARISON OF DYNAMICAL FEATURES 

In this section we go beyond the mean-field treatment of the previous sections to compare 
and contrast some of the dynamical features observed in long kinetic Monte Carlo simulations 
of Model A and Model B. For both models we performed multiple simulations of 2^^ = 
33 554432 generations j69|. with a genome of length L = 13 (2^'^ = 8192 potential species) 
and a mutation rate per individual of = 10~^. The fecundity F was set to 4 for Model A 
and 2 for Model B. 

The model parameters used in our simulations are chosen for computational feasibility 
with a view to keeping the system in the realistic regime where A/'(t) -C 2^ and A/'(t) <C 
Ntot{t) at all times. At the same time, we want to study the dynamics during the late-time, 
statistically stationary regime, which should therefore be reached in a reasonable time. As 
a result of these restrictions, L, R, and Nq are quite small, while fi is quite large, compared 
to natural populations. Results of tests for Model B with larger L and R and smaller 
are summarized in Appendix O Analogous, but less exhaustive tests for Model A (not 
reported in detail here, but see brief summary in Appendix [C]) yield similar results. These 
tests indicate that the products Rfi for Model B and Nq^i for Model A are proportional to 
the average number of mutant individuals produced per generation, A^tot/^- We also note 
that this number is half of the universal biodiversity number in Hubbell's neutral theory of 
biodiversity jlH]. Even though the species in the models studied here are strongly interacting, 
this parameter appears to play a similar role in determining the average diversity as it does 
in the neutral theory. Furthermore, the power laws that are observed and discussed in 



11 



Sec. IIV Bl are found to be robust, even against parameter changes that change 



A. Evolution of mean-field parameters in the transient regime 

The simulation runs were started with 100 individuals of a single species (or of a single 
producer species for Model B) and allowed to evolve under the dynamics described in Sec. [Tll 
At time points sampled at approximately constant intervals on a logarithmic scale, the 
"stable core" of the community was identified as follows. The fixed-point population sizes 
for the species present in the community at that time were calculated according to Eq. ([6]), 
and species corresponding to negative population sizes were excluded (starting from the one 
with the smallest simulated population size) until a stable, feasible community was obtained. 
This procedure was found to essentially correspond to removing species with populations 
below 100, which mostly were unsuccessful mutants of the "core" species. Once the core 
community was identified, the corresponding mean- field parameters B, £ (for Model B only), 
N^ot^ ^^cl $max were calculated from Eqs. ([8]) through f|TT]) . 

Typical time series of these mean-field parameters are shown in Figs. [3] and HI in which 
the early-time transient regime is emphasized by using a logarithmic time axis. For both 
models a clear tendency is seen for the mean-field parameters to increase during the transient 
regime, for so to fluctuate randomly about their long-time averages during the the later- 
time, statistically stationary regime. Both models thus evolve from an initial low-population 
community, toward a sequence of QSS communities that optimize their population sizes, 
subject to the constraints inherent in the respective models. 



B. Long-time dynamics in the statistically stationary regime 

Here we concentrate on comparing and contrasting the dynamical properties of the two 
models in the statistically stationary regime that follows the transient regime discussed 
above. The emphasis will be on the temporal fluctuations of the total population sizes and 
the diversities of the resulting communities. 

Typical time series of the total population sizes A^tot(^) are shown in Fig. [51 They are 
characterized by QSSs on different time scales, separated by periods of high evolutionary 
activity. In addition to the total population sizes, we also studied the diversities, defined as 
the number of major resident species. In order to obtain an approximation for the number of 
major, or "core" species (which can be thought of as the wildtypes in a quasi-species model 
[isl, ISOj), we filter out the low-population species that are most likely unsuccessful mutants 



of the wildtypes. This is achieved by using the exponential Shannon- Wiener diversity index 

D{t) = e^I{"^W}] , (15) 

where 

S[{njm = - E p/Winp/W (16) 

{l\pi{t)>0} 



with pi{t) = ni{t)/Ntot{t) is the information-theoretical entropy |54j . l55l |. Typical time 
series for D{t) in the two models are shown in Fig. [61 Just like in the time series for the 
total population sizes, the intermittent structure consisting of QSSs on different timescales, 
separated by periods of high evolutionary activity, is clearly seen. 
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FIG. 5: (Color online.) Typical time series of the total population sizes, Ntotit), for Monte Carlo 
runs of 2^^ = 33 554 432 generations each. To reduce the size of the figure file, the data are sampled 
only once every 2048 generations. While this slightly reduces the apparent range of the short-time 
fluctuations, the general shapes of the time series are preserved. The data are from the same runs 
shown in Figs. [3] and H (a) Model A. (b) Model B. 
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FIG. 6: (Color online.) Typical time series of the exponential Shannon- Wiener diversity index 
D{t) in the same 2^^-generations Monte Carlo runs shown in Fig. [5j The data were downsampled 
as in Fig. El (a) Model A. (b) Model B. 



Statistical information for characteristic time intervals that describe the dynamics can 
be extracted from our data. One such time interval is the duration of a QSS. One way 
to determine this is to use a cutoff on the magnitude of the diversity fluctuations, whose 
probability densities for the two models are shown in Fig. [7|(a). For both models, the 
probability densities consist of a Gaussian central part representing the fluctuations during 
the QSS periods {65-1 , flanked by "wings" that correspond to the large fluctuations during the 
evolutionarily active periods. We choose a cutoff i/c = 0.015 for Model A and 0.010 for Model 
B, in both cases corresponding to the transition region between the Gaussian central peak 
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FIG. 7: (Color online.) (a) Normalized histograms representing the probability density of the 
logarithmic derivative of the diversity, dS{t)/dt. The data were averaged over 16 generations in 
each run, and then averaged over 16 independent runs for Model A (solid, marked A) and 12 runs 
for Model B (dashed, marked B). The central parts of both histograms are well fitted by a Gaussian 
(dotted, marked G). (b) Log-log plot of normalized histograms representing the probability density 
of the durations of QSSs, estimated as the periods between times when |dS'(t)/dt| exceeds a cutoff 
Uc {dS{t)/dt was averaged over 16 generations as in part (a)). Model A (solid) and Model B 
(dashed). Results for Model A were averaged over 16 runs, and for Model B over 12 runs. The 
error bars are standard errors, based on the spread between runs. The two dot-dashed straight 
lines represent time~^ and time~^ power laws, respectively. 



and the wings. The duration of a single QSS then corresponds to the time interval between 
consecutive times when |dS'(t)/dt| exceeds i/c. Log-log plots of the resulting probability 
densities for the durations of QSSs in the two models are shown together in Fig. El^b). While 
both show approximate power-law behavior over five decades or more in time, there is an 
important difference: the exponent for Model A is near —2, while for Model B it is closer to 
-1. 

A different time interval of interest is the lifetime of a particular species, defined as the 
time between its origination and eventual extinction. Log-log plots of histograms of the 
species lifetimes in the two models are shown in Fig. [HI In contrast to the case of the QSS 
durations, the species-lifetime distributions are very close for the two models, both showing 
approximate time~^ behavior over near seven decades in time. The observed exponent is 
significantly different from —3/2, which would correspond to the simple hypothesis that the 
lifetime distributions simply correspond to the first-return-time distribution for a random 



walk of nj [40|, |41j. On the other hand, a lifetime distribution with an exponent of —2 
is consistent with a stochastic branching process [45j. (Simulations of neutral versions of 
both models also yield species lifetime distributions that decay approximately as time~^, 
but with a sharp cutoff near 10^ generations.) Lifetime distributions for marine genera that 
are compatible with a power law with an exponent in the range —1.5 to —2 have been 



obtained from the fossil record |41l . |42| . However, the possible power-law behavior in the 
fossil record is only observed over about one decade in time - between 10 and 100 million 
years - and other fitting functions, such as exponential decay, are also possible. Nevertheless, 
it is reasonable to conclude that the numerical results obtained from complex, interacting 
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FIG. 8: (Color online.) Log-log plot of normalized histograms representing the probability density 
of the lifetime of a particular species. Model A (solid) and Model B (dashed). Results for Model 
A were averaged over eight runs, and for Model B over 12 runs. The error bars were calculated as 
in Fig. W[i>)- The dot-dashed straight line represents a time~^ power law. 
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FIG. 9: (Color online.) Plots showing the labels I of the most highly populated species vs time. 
The data are for the same 2^^-generation Monte Carlo runs shown in Figs.[3H6l Gray (green online): 
nj G [101,1000]. Black: nj > 1001. (a) Model A. (b) Model B. For clarity, only producer species 
are shown for this model. The consumer species give a similar picture. See further discussion in 
the text. 

evolution models that extend over a large range of time scales support interpretations of the 
fossil lifetime evidence in terms of nontrivial power laws. 

The difference in the power laws for the QSS durations and the species lifetimes is a 
puzzling result. A likely explanation can be gleaned from the data shown in Figs. 191 and fTOl 
These show the species labels of highly populated species as functions of time for both mod- 
els. From Figs. [9](a) and [TOT a) it is seen that all or most of the horizontal lines representing 
populated species at a given time for Model A start and stop almost simultaneously, indi- 
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FIG. 10: (Color online.) Magnified version of Fig. [U showing times between 1.0 x 10'' and 1.6 x 10^ 
generations. In addition to the species shown in Fig. [U the ones with nj E [11, 100] are also shown 
(light gray, yellow online), (a) Model A. (b) Model B. See further discussion in the text. 



eating that species originations and extinctions in this model are highly synchronized. In 
other words: whole communities in Model A tend to go extinct and be replaced with an en- 
tirely new community within a short time. In contrast, in Figs. Mjo) and [TOT b) the different 
horizontal species lines for Model B stop and start at different times. This indicates that 
communities in this model are much more robust, and extinction events seldom wipe out 
more than a part of the total community. Thus, QSSs would be expected to be more long- 
lived (but also less clearly defined) for Model B, than for Model A. This theory is supported 
by the strikingly different QSS community structures produced by the two models, shown in 
Fig. [m The typical QSS community shown for Model A is a small cluster of mutualistically 
interacting species, while the typical community shown for Model B has the character of a 
simple food web. Extinctions in the latter are likely to be confined to a single branch of the 
web. A detailed discussion of the community structures generated by Model B, including 
comparisons with data from real food webs, is found in 47 . 

The arbitrariness inherent in the cutoff that must be used to extract QSS duration dis- 
tributions from fluctuations in the diversity or other time series can to some extent be 
eliminated by mapping distributions obtained with different cutoffs onto a common scaling 



function [43|, |49|. However, an analysis method that completely avoids cutoffs is that of 
calculating power spectral densities (the square of the temporal Fourier transform). Power 
spectra (PSDs) are therefore shown in Fig. [12] for both models. The PSDs for the diversity 
are shown in Fig. [T2](a). and for the total population size in Fig. [TW b). Although there are 
clear deviations, the overall behaviors for both quantities and for both models are compat- 
ible with a 1// power law over many decades in frequency. In the high-frequency regime 
the population-size PSDs have a significant background of noise, presumably caused by the 
rapid population fluctuations due to the birth and death of individual organisms. For very 
low frequencies there is little reason to believe that there should be large differences between 
the behaviors of the two quantities for the same model. We therefore think it is reasonable 
to consider the difference between the slopes of the diversity and population-size PSDs as 
an indication of the true uncertainty in the PSDs at the lowest frequencies. Better estimates 
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FIG. 11: (Color online.) Typical QSS core communities observed in simulations of the two models. 
The numbers are the species labels /, and the sizes of the circles represent the individual species' 
fixed-point populations, obtained from Eq. ([U]). A black arrow pointing toward species / indicates 
a positive Mjj, with magnitude indicated by the line thickness and size of the arrow head. The 
thin gray lines without arrowheads (red online) connect nearest neighbors in genotype space, (a) 
Model A. This particular QSS community is the one included in Table 2 of [(35] and line 10 of Table 
I of [sB]. (b) Model B. This is the food web representing the QSS community existing between 15 
and 18 million generations in the simulation shown in Figs. [H EJ^b), El^b), [9|^b), and 110(b). It has 
three trophic levels above the resource node. 

in this regime would require several orders of magnitude longer simulations. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have compared and contrasted the long-time evolutionary dynamics of 
two individual-based models of biological coevolution, using both analytic linear stability 
calculations and large-scale kinetic Monte Carlo simulations. These models involve universal 
competition and ignore important effects such as adaptive foraging. They are therefore 
not highly realistic. However, the fact that the results of numerical simulations can be 
compared with exact analytical results make these simplified models ideal as benchmarks 



for simulations of more realistic coevolution models in the future [48 . 

A central result of the analytic study is that, in the absence of mutations, the total pop- 
ulation size of a fixed-point community, N*^^, is described by a community fitness function 
that is cubic in A^tot- In addition to the present application, such a model is also applica- 
ble to nonequilibrium phase transitions in such diverse systems as epidemics, lasers, and 
autocatalytic chemical reactions. However, these evolution models differ from those kinds 
of systems by the important effect that, in the presence of mutations, the model parame- 
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FIG. 12: (Color online.) Log-log plots of PSDs of (a) the diversities D(t) and (b) the total 
population sizes Ntot{t) for Model A and Model B. The PSDs are averaged over each octave in 
frequency, and then over 16 runs for Model A and 12 runs for Model B. The error bars were 
calculated as in Fig. [TJ^b). The dashed straight lines represent 1/f power-law behavior. 



ters G and S are no longer externally imposed constraints, but rather evolve as far in the 
direction of positive G and large S as allowed by the internal constraints of the particular 
model. As a result. Model A, in which the elements of the interspecies interaction matrix 
M are randomly distributed on an interval that is symmetric about zero, evolves to produce 
communities that are heavily biased toward mutualism. The effective interaction variable 
G adjusts to a positive value. In contrast, the predator-prey Model B, in which the inter- 
species interactions are antisymmetric, is constrained to nonpositive values of G, for which 
a nonzero population size can only be sustained through the external resource R. These 
results are illustrated in Fig. [D^a). 

In a recent study of a similar, but different model for coevolution G^l, the emergence of 
strongly mutualistic communities from initially unbiased conditions was also observed. In 
that model, mutants are very similar to their parents, except for their interactions with a 
few other species ("local mutations"), and the authors suggest that the evolution of mutu- 
alism is related to this feature of their model. However, the mutations in our models would 
be "global" in the language of [iO], which leads us to the conclusion that the emergence of 
mutualism is common in models where direct mutualistic interactions are allowed. Rather 
remarkably, we have found little difference in the dynamics between the version of Model 
A studied in this paper and a version with strongly correlated interactions, and thus more 



"local" mutations [53|. It remains an important problem to reconcile this tendency for evo- 
lution of mutualism with the obvious requirement that biomass cannot be created without 
energy input. While predator-prey interactions are easy to reconcile with energy conserva- 
tion, direct mutualistic interactions are more difficult to interpret in an energy framework [61] 
(although effective mutualism is common in nature 0, 0, EI, Q). We believe this empha- 
sizes the importance of distinguishing between direct mutualistic interactions, as in Model 
A, and more realistic effective mutualisms, which merely mean that a pair of elements in 
the community matrix, Ajj and Ajj, are both positive. The unrealistic emergence of direct 
mutualism in Model A is shared by the Tangled-nature model 10, 11, 2^, of which it is a 
simplified version. 
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Beyond the mean-field studies and the simulation results for average population sizes, we 
have also studied the temporal fluctuations of both the diversity and total population size 
for Models A and B. We find that the probability distributions of the lifetimes of individual 
species in both models are very similar, showing power-law decay with an exponent near 
—2 over near seven decades in time, as seen in Fig. [SI This exponent value is consistent 
with some interpretations of the available data for the lifetimes of marine genera in the fossil 
record 41, 4^, but other interpretations of the fossil evidence are also possible. Similarly, 



power spectra for the diversity, as well as for the total population size, show reasonable 
(although not perfect) 1// behavior over many decades in frequency, as seen in Fig. [121 It 
is therefore very interesting that the probability distributions of the durations of individual 
QSS periods in the two models also both show reasonable power-law decay, but with different 
exponents: near —2 for Model A and close to —1 for Model B, as seen in Fig. [7|(b). This 
result, which we found quite surprising at first, makes sense in light of the observation that 
the extinctions of major species are highly synchronized in Model A, while they are much 
less so in Model B. While communities in Model A tend to collapse completely when an 
aggressive mutant arrives and/or a major species goes extinct, communities in Model B are 
much more resilient and extinctions most often only extend to one or a few branches of the 
resident food web. This effect is illustrated in Figs. [9], [TUl and[TTJ The long-time correlations 
that give rise to these extended power laws are clearly due to the interspecies interactions. 
Indeed, our exploratory simulations of neutral versions of both models (not shown here) 
exhibit no correlations beyond about 10^ generations. 

Our observation of the high resilience of Model B against complete extinction of communi- 
ties is consistent with observations of extinction avalanches of limited size in the Web-world 
model by Drossel et al. lB|, who argue that their model is therefore not self-organized crit- 



ical. Together with our observation of the self-optimization of the evolution models studied 
here to points away from the transcritical bifurcation point, these observations may support 
a conclusion that models of coevolution that take reasonable account of the dynamics at 
the ecological level (even if they are extremely simplified) are not in general self-organized 
critical. Such a conclusion, which appears reasonable, would be in disagreement with a 
number of recent theories of extinction (sl. 41. 43l . 

The discussion in the preceding paragraphs brings out some of the differences and simi- 
larities between the models studied here, and some of the other physics-inspired evolution 



models that have been introduced in recent years. Like the Tangled-nature model [10|, [III, 129 
our models are individual-based, and Model A also shares with that model the unrealistic 
emergence of mutualistic communities. The power-law behaviors are also reminiscent of 
species-based (rather than individual-based) extremal-dynamics evolution models like the 
Bak-Sneppen model and its generalizations |3l. [41. Si, and in [4§| we speculated that Model 
B might correspond to a zero-dimensional extremal-dynamics model fi^- However, more 
accurate exponent estimates for Model B 47(1 , as well as the above comparison of Model B 



with the Web-world model [l6|, make this conjecture less likely to hold. 

In conclusion, the results presented here indicate that further work on models of 
macroevolution that are based on events on ecological time scales, with comparisons of 
the results with data from the fossil record, as well as from laboratory experiments and 
extant food webs, is highly desirable. In a separate paper we consider in detail the structure 
and dynamics of the food webs that develop within Model B, and we compare these with 
data for real food webs [4j 



19 



Acknowledgments 

The author thanks R. K. P. Zia, B. Schmittmann, P. Beerh, G. J. P. Naylor, and V. Sevim 
for useful discussions and comments on the manuscript, and V. Sevim for the data on 
Model A included in Fig. [HI He also gratefully acknowledges hospitality at Kyoto University 
by H. Tomita in the Department of Fundamental Sciences, Faculty of Integrated Human 
Studies, and H. Fujisaka in the Department of Applied Analysis and Complex Dynamical 
Systems, Graduate School of Informatics. 

This work was supported in part by U.S. National Science Foundation Grants No. DMR- 
0240078 and DMR-0444051, and by Florida State University through the School of Com- 
putational Science, the Center for Materials Research and Technology, the National High 
Magnetic Field Laboratory, and a COFRS summer-salary grant. 



APPENDIX A: DERIVATION OF THE COMMUNITY FITNESS FUNCTION 



In this Appendix we provide a conventional derivation of the cubic form of the community 
fitness function $(A^tot) in a simple mean-field approximation. The derivation provides an 
explanation for the prefactor (l — l/F) in Eq. (ITT]) , as well as intuitively clear approximations 
for the coefficients G and S. It also provides justification that the equation to be integrated 
to obtain $(A^tot) is indeed Eq. ([7]), rather than this equation multiplied or divided by some 
power of A^tot- The derivation is based on the time-dependent Ginzburg-Landau equation for 
a system with nonconserved order parameter 27|, |3l| , which for our current systems takes 
the form 

dt -dN,^,- ^^^^ 

Identifying dni{t)/dt with nj{t + 1) — njit), we obtain from Eqs. ([IHS]) in the absence of 
mutations: 



dni{t) 
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Expanding this nonlinear equation of motion around its fixed point, we get 
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To obtain the simplest mean-field approximation for dN^ot/dt (exact for A/" = 1), we set 

ni ^ N,ot/M, bj ^ M~^Y!ibi = (b), VI - ^-'E'lVi = iv), and Mu ^ Af-^E'ijMij = 
(M) , where the primes on the sums indicate that they are restricted to the Af species with 
rii > 0. This yields 
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Integrating the right-hand side as prescribed by Eq. (lAll) . we find 
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This result has the same cubic form as Eq. (fTTI) . with the approximate coefficients {if) S 
and ((M) — (6)) ^ 6. The exact fixed-point solution for iVtot requires the use of £ and 6, 
but the approximate forms obtained here may provide a more intuitive understanding of the 
significance of these coefficients. 



APPENDIX B: DERIVATION OF OPTIMUM POINT FOR MODEL A 



We can obtain a theoretical estimate for the point (O, A^tot) in Model A. We assume for 
simplicity that all the off-diagonal Mjj have the same value, a. Then, using the definition 
of G, and remembering that F = A for Model A, so that 6/ = — In 3 for all /, one can show 
that = (1 — l/A/')a + ln3, where M is the total number of species in the community. This 
yields the absolute maximum value for B equal to (1 -|- InS) 2.10 for M = oo (shown by 
vertical full lines in Fig. [1]). However, it has been shown [sO] that the most probable number 
of species in a community within this model is finite and given by 

A/"^ ^Lln2/ln(l/g) , (Bl) 

where q is the probability of finding a pair of interactions, M/j and Mj/, conducive to this 
community. With the Mjj independently distributed on [—1, +1], the probability of drawing 
a pair that are both larger than some given value m G [— 1,+1], is g = [(1 — m)/2]^. In 
our approximate formula for 0, we now replace M by N"^ with this value of g, and a by the 
average of a variable uniformly distributed over [m, 1], which is (1 + m)/2. This yields 
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+ ln3 (B2) 



for L = 13. This is a concave function of m with a single maximum at mmax, which can be 
found numerically. The result is mmax ~ 0.512, which yields 6max ~ 1.618, A'tot ~ 3236, and 
Mjj = (1 -|- minax)/2 ~ 0.756. These values are in excellent agreement with those obtained 
from the Monte Carlo simulations, and they are shown as gray dots (turquoise online) in 

Fig.m 



APPENDIX C: TRIALS WITH OTHER PARAMETER VALUES 

In addition to the parameter set used in the main simulations presented above, we also 
performed trial simulations for Model B with L = 20 (1048 576 potential species), /i = 
2.5 X 10~^, and R = 8000. As each simulation run with R = 2000 took about two weeks of 
CPU time, and each run with R = 8000 took at least four weeks, relatively few trial runs 
were performed. No qualitative differences were observed in the quantities reported on the 
basis of the main simulation series. 

To obtain estimates of the effects of L, fi, and R on the the mean stationary levels of the 
total population size, A^tot, and diversity, D, and the mean time to reach stationarity, r, we 
fitted these variables by the exponential function a[l — exp(— t/r)]. Numerical results for r, 
D, and A'tot, based on twelve runs for L = 13, R = 2000, and fi = 10~^ (the original runs 
used in the main part of this paper); twelve runs for L = 20, R = 2000, and fi = 10~^; nine 
runs for L = 20, R = 2000, and = 2.5 x 10"^ and five runs for L = 20, R = 8000, and 
jj, = 2.5 X 10 ~^ are compiled in Table [B 
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TABLE I: Numerical results for the main simulations of Model B with L = 13, compared with trial 
runs with L = 20 and varying R and /i. The unit for R and A'^tot is individuals, for ^ it is mutations 
per individual offspring per generation, for r it is 10^ generations, and for D it is species. Thus, 
Rfx is proportional to the total number of mutated offspring per generation. The uncertainties are 
standard errors, based on the spread between individual runs. cJr is the standard deviation of r 
over the independent runs. A ratio t/iTt- 1 indicates that r may be approximately exponentially 
distributed, while the even lower ratio for = 0.5 indicates very large variations from run to run. 



L 
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Ai 


R^l 


Runs 
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t/(Jt 


D 


Atot 


13 


2000 


10-3 


2.0 


12 


0.35±0.07 


1.08 


10.9±0.3 


7700±300 


20 


2000 


10-3 


2.0 


12 


0.8±0.2 


0.93 


14.7±0.3 


13 700±400 


20 


8000 


2.5 X 10-^ 


2.0 


5 


1.2±0.3 


0.90 


20.4±0.6 


58 000±7000 


20 


2000 


2.5 X 10-^ 


0.5 


9 


5±2 


0.53 


10.3±0.2 


12 000±2 000 



1. Effects of the genome size 

With L = 13, the whole pool of potential species is visited in a typical 2^^-generation 
simulation, while only about 10% are typically visited for L = 20 with /i = 10-3 and 
R = 2000. However, the "revivals" of extinct species, seen frequently with L = 13 [sO] and 
much less so with L = 20, has no effect on the observed power laws. Twelve full runs were 
performed for this parameter set. All the subsequent trial simulations were performed with 
L = 20. _ 

As expected from Eq. fIBip . D appears to be proportional to L. The significant increase 
in Atot when L is increased from 13 to 20 (not seen in trial runs with L = 20 for Model A) is 
probably due to the increased probability of finding species with smaller bj and thus closer 
to the population divergence for Model B at B = (see Fig. [T](a)). 

2. Effects of reduced jj. and increased R 

We found r to be approximately inversely proportional to fi and R, while the mean di- 
versity increases (but apparently sublinearly) with both /i and R. Within the uncertainty, 
r was the same for the diversity and population size, and the reported results are there- 
fore based on both quantities. The mean total population size Atot appears to be roughly 
independent of n and approximately linearly dependent on R. This latter proportionality 
(which is also expected from the analytical result, Eq. ([9])) means that Rfi oc Atot A* , the av- 
erage total number of mutant individuals per generation. Even though the models studied 
here have strong interspecies interactions, we note that Atot/^ is one half of the fundamental 
biodiversity number in Hubbell's neutral theory of biodiversity [l^]. This parameter appears 
to play an analogous role in determining the diversity in the present models as it does in 
the neutral theory. 

However, more important than the above results is that none of the observed power 
laws changed under parameter variation within this range. In particular, the exponent for 
the PDF of individual species lifetimes remained near —2 for all parameter sets, while the 
exponent for the QSS duration PDF remained close to —1. In fact, allowing for the larger 
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overall intensity in the population PSD for R = 8000 and to within the numerical accuracy, 
both PDFs and PSDs for all three parameter sets with L = 20 can be overlaid graphically 
with the ones presented for L = 13, fj, = 10~^, and R = 2000 elsewhere in this paper. Based 
on these trials we therefore conclude that the parameters used in the main study are also 
representative for systems with larger L and R and smaller /i and are well chosen to study 
the stationary dynamics of the models on long time scales. 

Analogous, but less exhaustive, tests for Model A give similar results. Runs with Nq = 
16 000 and fi = 1.25 x 10~^ {Nqh = 2) gave A'tot ~ 1.5 Nq as in the main simulations with 
A^o = 2000, and D ^ 3.6, a statistically insignificant increase from the value of approximately 
3.3 for the simulations with A'o = 2000 and fi = 10~^ 13 ■ 
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